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In  this  report  we  describe  progress  in  developing  and  implementing  a 
general  theoretical  approach  to  describing  the  properties  of  defects  and 
Impurities  of  a  general  nature  in  non-metallic  solid  systems.  This  approach 
combines  fully  correlated,  fully  self-consistent  electronic  structure 
determination  of  the  electrical  and  mechanical  properties  associated  with 
neutral  or  charged  defects/impurities  in  or  on  a  non-metal.  The  system 
remote  from  the  defect  is  described  by  the  shell  model  which  incorporates 
self-consistentlv.  host  polarization  and  distortion.  This  results  in  our 
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being  able  to  obtain  absolute  energies  of  the  impurity  ions  in  the 
host  and  their  interaction.  The  model  is  free  of  adjustable  or  undefined 
parameters.  This  project  is  of  non-trivial  magnitude  and  the  current 
computer  implementation,  which  is  functional  in  our  laboratory,  consists 
of  a  program,  ICECAP,  which  is  about  100,000  statements  long.  This  pro¬ 
gram  is  the  result  of  extensive  collaboration  between  our  group  and  that 
of  Prof.  J.  M.  Vail,  University  of  Manitoba,  and  of  Dr.  A.  M.  Stoneham, 
Harwell,  AERE. 
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In  this  report  we  describe  progress  in  developing  and 


implementing  a  general  theoretical  approach  to  describing  the 


properties  of  defects  and  impurities  of  a  general  nature  in  non- 
metallic  solid  systems.  This  approach  combines  fully  correlated, 
fully  self-consistent  electronic  structure  determination  of  the 


electrical  and  mechanical  properties  associated  with  neutral  or 
charged  defects/impurities  in  or  on  a  non-metal.  The  system 


remote  from  the  defect  is  described  by  the  shell  model  which 


incorporates  self-consistently,  host  polarization  and  distortion. 
This  results  in  our  being  able  to  obtain  absolute  energies  of  the 
impurity  ions  in  the  host  and  their  interaction.  The  model  is 


free  of  adjustable  or  undefined  parameters.  This  project  is  of 
non-trivial  magnitude  and  the  current  computer  implementation, 
which  is  functional  in  our  laboratory,  consists  of  a  program, 
ICECAP,  which  is  about  100,000  statements  long.  This  program  is 
the  result  of  extensive  collaboration  between  our  group  and  that 
of  Professor  J.  M.  Vail,  University  of  Manitoba,  and  of  Dr.  A.  M. 
Stoneham,  Harwell,  AERE. 
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The  aim  o-f  this  research  line  is  to  obtain  absolute  energies 
of  impurity  ions  or  defects,  including  charged  states,  jn  non- 
metals.  In  addition,  the  absolute  excitation  energies  of  low 
lying  excited  states  are  studied.  These  calculations  will  not 
simply  be  of  one-electron  energies,  nor  will  spectroscopy  be 
obtained  by  using  Kaopmans’  theorem,  but  rather  by  evaluation  of 
total  system  energies  and  their  differences.  These  calculations 
will  include  electron  correlation  corrections  and  multiplet 
splittings.  Due  to  using  differences  in  total  energy  to  generate 
spectroscopic  information,  coulomb-hole  attractions  will  be 
included  directly  and  not  as  a  perturbation.  i 

Such  studies  are  important  in  a  wide  range  of  applications, 
such  as  i)  spectroscopic  and  laser  applications,  where  it  is 
necessary  to  know  about  the  stability  of  a  promising  charge  state 
of  a  particular  impurity  in  a  novel  host,  ii)  solid-state 
reactions,  where  gas  sensors  and  oxide  reactions  serve  as 
examples,  in  which  the  impurity  acts  as  a  source  or  sink  for 
electrons,  exchanging  carriers  with  other  species  at  the  surface 
or  in  the  bulk,  and  iii)  stability  of  charged  imperfections  in  or 
on  unstable  solids  (energetic  solids),  where  local  deviations  from 
periodicity  may  control  stability  of  entire  systems,  iv)  etc. 

It  is  clear  that  these  examples  show  only  a  small  fraction  of 


the  variety  of  important  properties  which  depend  crucially  upon 
the  absolute  energy.  We  observe  that  calculations  of  the  sort 
described  lead  naturally  to  prediction  of  a  wide  variety  of 
related  quantities.  Charge-state  stability,  for  example,  is  one 
aspect  of  the  prediction  of  photoionization  energies,  but  also  can 
lead  to  the  identification  of  which  of  several  alternate  sites  in 
a  complex  substance  is  the  stable  one,  to  the  identification  of 
charge-compensation  mechanisms,  and  to  the  prediction  of  limiting 
solubilities  with  greater  accuracy  than  current  crude  charge  and 
size  misfit  models. 

A  second  key  output  of  such  studies  can  be  the  development  of 
interatomic  potentials,  which  have  a  sound  theoretical  basis, 

»  <p 

using  the  best  available  current  technique.  These  calculations 
must  therefore  include  electron  correlation  explicitly.  This 
aspect  shall  emphasize  those  cases  for  which  empirical  potentials 
are  problematic.  These  systems  include  systems  as  0*-,  S*-,  Se*-, 
Te*-,  etc.  In  these  cases  the  free  double  negative  ion  is  not 
stable  in  free  space  but  is  only  stabalized  by  its  environment. 

The  ionic  polarizability  of  these  entities  is  highly  dependent 
upon  the  host  and  perhaps  even  the  local  site  occupied  in  the 

I 

host. 

A  second  category  includes  systems  for  which  empirical  forms 

are  not  available.  This  may  be  because  the  species  do  not  occur 
at  the  proper  separation  in  wel 1 -documented  perfect  crystals,  or 
because  special  species  are  involved,  such  as  F-  -Fa-  or  I~  -0*-, 

or  because  unusual  cation  -  cation  interactions  occur.  The  last 

case  is  problematic  since  cations  are  rarely  found  close  enough  to 

each  other  to  yield  more  than  a  very  small  component  of  the  total 


energy  of  perfect  crystals. 

A  third  category  includes  a  wide  variety  of  other  cases  which 
include  systems  in  unusual  charge  states  such  as  FeB+  in  SrTiO»  or 

Fe+  in  MgO.  Impurity  systems  in  their  excited  states  are  also 
included  here,  which  are  not  included  in  traditional  approaches. 

Any  interatomic  potentials  developed  by  this  project  are 
subject  to  verification  in  three  ways.  These  include:  i)  checks 
that  these  potentials  provide  good  values  in  applications  to 
defect  properties  (thermodynamic  vibrational,  lattice  expansion, 
etc.)  parallel  to  similar  studies  produced  previously,  ii)  checks 
of  predictions  over  a  range  of  crystals  for  properties  of  wider 
interest,  such  as  electro  optic  constants  in  complex  oxides,  to 
see  if  trends  are  accurately  reproduced,  iii)  checks  of 
transferability  for  species  such  as  0*-  and  Mg5**-  to  see  to  what 
extent  their  interatomic  potentials  are  the  same  in  MgO  and  in 
complex  oxides  such  as  BMAG  (Ba*  Mg  Gea  0*-)  or  in  the  bulk  and 
near  a  surface.  This  latter  study  is  particularly  interesting  for 
the  double  negative  ions  in  that  they  are  not  stable  in  free 
space. 

These  studies  are  directly  applicable  to  a  wide  variety  of 
energetic  solids  including  (but  not  limited  to)  CH^,  (NO)^* 
Nitromethane,  RDX. 


It  is  of  practical  importance  to  be  able  to  calculate  reliably 
properties  of  defects  in  crystalline  materials.  This  work  relates 
to  methods  developed  for  point  defects  in  ionic  crystals,  for 
which  many  successful  ground  state  calculations  have  been 
performed  on  the  basis  of  a  classical,  discrete-ion  model: 


dipole  polarizable  point  charges.  Each  ion  is  represented  as  a 
point-charge  care  coupled  harmonically  (force  constant  K)  to  a 
uniformly  charged  (charge  Y>  massless  spherical  shell  of 
indeterminate  radius,  and  ions  interact  through  Coulomb  and  short- 
range  shell-model  potentials,  v(r>.  The  latter  are  exemplified  by 
the  Buckingham  forms 

v(r)  =*  B  exp(-r/p)  -  (C/r*).  (1) 

The  use  of  this  form  is  not  essential. 

In  the  defect  cluster,  excess  electrons  are  treated  quantum- 
mechanical  ly, .usually  incorporated  with  the  electrons  of  some  near 
neighbour  ions.  In  our  work,  the  electrons  of  the  cluster  are 
treated  in  the  unrestricted  Hartree-Fock  (UHF)  self-consistent 
field  (SCF)  approximation  C31,  based  on  linear  combination  of 
atomic  orbital  (LCAO)  molecular  orbitals  (MO)  corrected  for 
correlative  effects  by  use  of  Many  Body  Perturbation  Theory 
(MBPT) .  In  most  non  metals,  the  electronic  structure  of  the 
ions/atoms/molecules  is  assumed  to  be  wel 1 -localized  about  the 
nuclei.  However,  this  will  not  automatically  follow  for  doubly 
negative  ions  on  the  outer  boundary  of  a  UHF-SCF  region  unless  the 
spatial  range  of  the  basis  AOs  is  restricted:  the  Coulomb  field 
of  the  surrounding  point-charge  shell -model  lattice  does  not 
impose  localization  upon  the  quantum-mechanical  region.  This 
difficulty  is  currently  overcome  by  associating  completely  frozen 
shell  ions  or  by  complete-ion  pseudopotentials  with  ions  that 
surround  the  UHF-SCF  region,  retaining  their  classical  dipole 
moments  and  shell-core  interaction  energies.  A  UHF-SCF  program 
C41  which  can  incorporate  such  pssudopotentials  (as  complete-ions 


or  as  ion  tores)  is  applied  in  this  work.  More  recent  theoretical 

developments,  not  yet  implemented  in  computer  code,  solve  the 

boundary  layer  problem  in  a  quantum  mechanically  exj|ct  way  as 

detailed  at  the  end  of  this  section. 

The  defect  cluster  embedded  in  a  shell -model  lattice  is 

described  mathematically  in  terms  o-f  lattice  (ionic  shell  and 
core)  coordinates,  collectively  denoted  R,  referred  to  as  the 

lattice  configuration,  linear  coefficients  in  the  LCAO-MO 
formulation  of  the  UHF-SCF  approximation,  collectively  denoted  £, 


referred  to  as  the  electronic  configuration,  nuclear  coordinates 


(nuclear,  pseudapotential ,  and  shell-model)  in  the  defect  cluster, 
collectively  denoted  gc,  referred  to  as  the  cluster  configuration 

and  coefficient  A,  of  excited  determinants  found  by  MBPT.  (The 


cluster,  and  therefore  Rc,  may  include  shell -model  ions  whose 


positions  are  anharmonical ly  perturbed  by  the  defect).  The  total 
energy  of  the  defect  E(Rc,g,g)  is  minimized  with  respect  to  R^c, 
and  JR: 


SE 
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yielding  a  variational  estimate  of  cluster  (Rc)»  electronic  (c>, 
and  lattice  (JR)  configurations,  and  of  total  energy  E  and  the 
electronic  wave  function  y(£). 

These  two  approaches  are  merged  in  such  a  way  that  in 

practice,  the  excess  electrons  and  UHF  and  pseudopotential  ions 

•re  first  simulated  by  fixed  point  charges,  from  which  HADES 

determines  the  polarized,  distorted  lattice  configuration  R,  and 

the  total  energy  E  of  the  shell-model  lattice.  The  shell-model 
point  charges  of  the  Lattice  (R) ,  along  with  the  nuclei  and 

pseudopotentials  of  a  fixed  cluster  configuration  Rc,are  now 
applied  as  a  background  potential  for  the  UHF-SCF  solution 


c,  yielding  total  electronic  energy  and  wave  -function  y(£>. 
Ideally,  the  point-charge  simulation  of  the  cluster  in  the  HADES 
calculation  should  have  all  its  electric  multipole  moments 
identical  to  those  of  the  UHF-SCF  cluster,  but  this  is  obviously 
not  practical.  One  therefore  matches  only  a  finite  set  of  low- 
order  multipoles.  This  is  accomplished  by  introducting  additional 
point-charge  simulators  into  the  HADES  calculation,  representing  a 
small  dipole,  quadrupole,  octupole,  etc.,  correcting  for  the 
discrepancies  between  HADES  and  UHF-SCF  up  to  a  given  multi pole 

order,  and  then  iterating  the  HADES/UHF-SCF  up  to  a  given 
consistency.  This  procedure  is  repeated  for  each  increment 

the  cluster  configuration,  until  minimum  total  energy  E  is 

obtained.  This  energy  iss 

6  -  <EH  -  E’c  *  Ei>  ♦  <Efl  +  <Efl  ♦  Ec  *  Ed*’  ,3, 

where  E,,  and  E.  are  defined  above,  E’  and  E’  are  the  shell -model 
HA  ’  c  s 

Coulomb  and  short-range  interactions  of  the  cluster  region  in  the 

HADES  calculation  that  will  be  replaced  in  the  UHF-SCF  calculation 

by  electronic  Coulomb  and  short-range  interactions  in  E^  and  by 

nuclear  and  pseudopotential  Coulomb  interactions  Ec,  and  E^ 

corrects  for  the  energy  of  the  dipoles  from  the  HADES  calculation 

that  become  associated  with  complete-ion  psuedopotentials. 

In  preliminary  work,  the  Lattice  R  includes  only  ions  out  to 
a  distance  beyond  which  the  defect  cluster’s  electric  monopole 

moment  (total  charge)  dominates,  where  HADES  uses  continuum  theory 

to  determine  the  discrete-lattice  polarization.  Since  the 

potential  of  such  a  polarized  continuum  would  be  constant  within 

the  cluster,  as  assumed  in  our  calculation,  there  is  a  small 

discrepancy  because  the  medium  is  in  fact  discrete.  This  will  be 


el iminated • as  time  permits. 

The  UHF-SCF  program  C43  used  in  our  work  has  available  a  many- 
body  perturbation  theory  correlation  correction  C53.  Because  UHF- 
SCF  calculations  become  very  time-consuming  as  the  number  of 
electrons  rises  above  200,  the  cluster  size  is  somewhat 
restricted,  and  correlation  correction,  which  is  additionally 
time-consuming,  cannot  be  applied  with  abandon.  Presently, 
clusters  of  more  than  a  few  lattice  spacings  in  radius  cannot  be 
analysed  in  this  way,  and  therefore  diffuse  electronic  states  of 
localized  defects  C63  are  not  accessible. 

The  physical  model  described  above  has  a  quantum-mechanical 
defect  region  (UHF-SCF  with  pseudapotentials  and  correlation 
correction),  with  perfect  lattice  boundary  conditions  (complete- 
ion  pseudopotentials  with  dipole  correction) ,  embedded  in  a 
classical  lattice  (shel 1 -model ) ,  solved  variationally  by  energy 
minimization  with  respect  to  cluster,  electronic,  and  lattice 
configurations,  to  finite-order  multipole  consistency.  CFig.  13 

The  computer  program  is  based  on  HADES  C13  and  UHF  C43 
programs,  both  of  which  have  been  extensively  tested,  refined,  and 
applied.  Consequently  any  point  defect  conf iguration  in  any  ionic 
crystal  host  lattice  geometry  can  be  analysed,  provided  shell- 
model  parameters  and  adequate  computer  time  are  available.  Atomic 
orbital  sets  may  include  s,  p,  d,  and  f  types  (except  that  the 
pseudopotential  option  cannot  accept  f-type  orbitals  yet).  Either 
norm-conserving  BHS  C73  or  Phi  1 1 ips-Kl einman  C83  pseudopotentials 
can  be  used,  and  correlation  corrections  can  be  included  in  the 
energy.  Octupole  consistency  is  presently  available.  The  program 
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is  currently  operating  on  the  Harwell  IBM  3081  computer,  and  the 
VAX-FPS  system  at  Michigan  Technological  University. 

The  dominant  computer  time  bottleneck  is  in  the  UHF-SCF 
process  end  in  application  of  MBPT.  Me  consider  now  some  computer 
dependent  ramifications  of  this  and  indicate  approaches  to 
solution. 

In  all  our  studies  the  Unrestricted  Hartree  Fock  (UHF)  method 
is  employed  as  a  starting  point,  as  is  the  normal  non-relati vistic 
Hamiltonian.  It  is  assumed  the  nuclei  are  infinitely  massive  and 
the  Born-Oppenheimer  approximation  is  used.  Ideally,  one  would 
like  to  solve  the  n-electron  Schrodinger  equation.  However,  exact 
solutions  to  (2)  are  seen  as  impractical  for  these  systems,  and  we 

resort  to  the  UHF  approximation.  That  is,  one  approximates  the 

solution  by  a  single  Slater  determinant  of  one  electron  orbitals, 

a 

*  i  In  the  UHF  approximation,  these  orbitals  are  constrained  to 
form  an  orthonormal  set  and  to  be  eigenstates  of  the  z  component 

of  spin.  The  orbitals  are  not  constrained  to  be  double  occupied 

or  to  have  well  defined  symmetry  properties.  Choosing  the 

orbitals  variational ly  yields  the  Hartree-Fock  equations 

F(P°>  »  €®  *?, 

1  i  x  ’  t 

(4) 
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P<x*,h>  is  the  operator  which  replaces  coordinate  x  with  x*. 

For  a  solid  system  with  low  symmetry,  such  as  a  solid  with  a 
point  defect,  solutions  to  even  the  UHF  system  of  equations  are 
expensive  to  encompass.  It  is  useful  to  make  use  of  the 
arbi trariness  of  the  Fock  equation  and  to  rotate  to  local 
solutions  if  possible.  The  way  this  is  done  has  been  given  by 
Kunz  and  Klein  C33  and  further  developed  by  Kunz  C53.  This 
technique  is  most  useful  for  non-metallic  systems  such  as  are 
studied  here.  One  formally  partitions  the  system  into  two  parts, 
the  cluster  to  be  studied  and  its  environment,  which  is  found 
using  the  HADES  method.  The  cluster  is  solved  self -consistently 

i 

in  the  field  of  the  environment. 

The  UHF  method  omits  correlation  effects.  A  brief 
.description  of  the  methods  being  currently  employed  is  in  order. 
Correlation  methods  to  be  used  for  extended  systems  are 
constrained  by  size  consistency  consi derati ons  E9, 103.  Our  group 
has  chosen  to  use  those  based  upon  multi -reference  many  body 
perturbation  theory  (MR-MBPT) .  Let  the  exact  Hamiltonian  be 
partitioned  into  a  "simple  Hamiltonian,"  Hq,  chosen  to  be  the  sum 
of  the  one-body  Fock  operators  for  the  n-body  system.  Thus 


t  . 


H*  H  +  V 
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and 
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Consider  the  -first  n  eigenstate  oif  Hq  separately.  There  may  be  no 
state  of  Hq  degenerate  with  these  n-states  unless  it  is  also 
included  in  the  n.  P  is  a  projector  onto  the  space  of  these  n 
states  and  is, 

P  *  £  UiX4i  1 . 
i=l  1  1 

Consider, 

<H  +  V)  ip 
o 

assume  we  wish  to  find  state  say.  Then 


H  *  »  (E  -  V)  ♦  and 
o 


(1  -  P)  (H  -  w,  )  y  =  <1->P>  (E-w.  -  V) 

o  l  l 

Commuting  (1  -  P)  with  <Hq  -  w^  permits  one  to  obtain  the  formal 
solution  for  ^ 

P*  -  ♦  -  <Hq  -  Wj)"1  (1  -  P)  (E  -  Wj  -  V) 

Now 

n  .  n 

P*  =  E  4  <4  U>  =  E  ir.4  .  =  4. 

j=l  J  J  j*l  J  J 


4  is  of  course  unknown.  Nevertheless  one  finds 


where 


T  =  <  1  -  <H  -  w,)-1  (l-P)  (E  -  w,  -  V)}-1' 
o  l  l 


One  may  obtain  the  energy  -from  the  secular  equation: 


<9> 


n 


V. 


(E  -  W.  >TT.  =  E  TT.  ik, 
11  k=l  k 


(10) 


where 

(11) 

From  a  utilitarian  point  of  view  these  equations,  (8)  -(11), 
are  not  -final  in  that  the  unknown  energy  E  occurs  in  the 
denominator  of  Eq.  (9).,  and  results  in  size  consistency  problems 
unless  treated  properly.  This  difficulty  may  be  circumvented 
here,  as  in  Rayleigh-Schrodinger  perturbation  theory,  by  using  the 
first  order  approximation  to  the  energy.  To  do  this,  and  to  solve 
these  equations,  one  expands  the  inverse  in  Eq.  (9)  in  a  power 
series 

T  -  1  +  (H  -  W, )_1  (1  -  P)  (E  -  W,  -  V)3  +  - 

o  1  1 


V. 
i  J 


<♦ 


VT,*j> 


+  C(H  -  WJ  1  (1  -  P)  (E  -  W,  -  V)]*^  -  (12) 

o  l  l 

In  this  case  the  first  approximation  to  E  is  found  by  solving 
<E  -  Wj)wt  -  ^  irk  <*,|v|*k>. 

(13) 

If  this  prescription  is  followed,  if  n  =  1  for  example,  one  simply 
recovers  ordinary  Rayleigh-Schrodinger  Perturbation  theory. 


Further  ramifications  of  the  use  of  MR-MBPT  and  ordinary  MBPT 


have  also  been  given  by  the  group  of  Bartlett  Cll,  123 


In  this,  as  in  most  other  numerical  studies,  the  algorithms 
chosen  are  designed  first  to  achieve  a  desired  level  of  precision 
and  only  then  chosen  to  maximize  efficiency.  After  all  there  is 
little  value  in  achieving  incorrect  results,  no  matter  how 
quickly.  In  these  studies,  we  fallow  one  of  the  conventional 
wisdoms  of  quantum  chemistry  and  expand  our  orbitals  in  a  basis 
set  of  gaussian  orbitals.  The  primitive  gaussian  orbital  is  of 
the  forms 

Xlx  i  j  k  ^  )  ■  <x  y  i  i+i+k)exp  (a^r-i^)2) 

( 14> 

This  function  has  the  advantage  that  all  necessary  integrals  over 
these  basis  functions  can  be  evaluated  in  closed  form  C103.  This 
allows  one  to  know  all  integrals  to  arbitrary  precision.  Using 
this  set,  one  may  construct  a  "contracted"  set  of  basis  functions, 
j  k ,  where 


tl  i  J  k  (r  -  f$»)  =  E  Aol  *1  <xi  j  k  (r  -  &  )  . 

1 


a 


<  IS) 


The  A’s  in  equation  (15)  are  assumed  given.  One  expands  the  Fock 

* 

solutions  in  terms  of  these  functions, 

,mP , 


*«  -  *i,  j,  k,  i  crj  k  1  «i  t  j  k 


(16) 


The  coefficients,  Ci  j  k  1 , .  are  found  using  the  Roothaan  method, 
which  is  in  reality  a  simple  linear  variation  C113.  It  is  this  we 

wish  to  discuss.  The  variation  is  performed  by  recalling  for  a 

particular  iteration,  the  say. 


P  (x  x1)  =  E  *  <x>  ♦  ^  <x  •  > 

m  m 

m=occ 


3  iE- 


...  5-  ,  .  (,  .  .  .  (x)  i  .  ( .  x  •.  > ,  E 
jkl  1  •  i  •  j  •  k  •  m 


r  mP  _  mP* 
ui  jkl  •  j  •  k » 1  • . 

The  sum  over  m  is  restricted  to  occupied  orbitals  and  mapping 
indices  i  j  k  1  into  an  index  i,  and  so  -forth,  so  that 

P(x  x‘>  -  E  f  .  <x)  *  ,.<x‘)  f 

i  i  »  i  • 


(17) 


where 


cP  _  _  -mp  _mP* 

^ii.  '  l  Ci  Ci. 

Ill 


P 

Thus  only  S^*  changes  from  iteration  to  iteration. 

The  Fock  problem  then  reduces  to  a  matrix  problem  (for  each 
iteration)  of  the  form 


F  ♦  =  E  D  ♦ . 


(18) 


The  matrices  D  and  F  are  given  as 


I 

l{  j> 
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Fij  =  <a|F|fJ> 
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(20) 


Clearly  the  needed  integrals  merely  need  be  evaluated  once  as  they 
don't  change  -from  iteration  to  iteration.  The  only  iteration 
dependent  quantity  is  Skl.  Let  us  define; 

I  h*  ea  I 

f  »  <i  I - Vs*  -  2  - - ||i> 

and  iJ  11  2m  il?-Ril' 


.Then 


9ijkl 


e* 


Fij  *  *ij  +  kl  l_9ijkl  9ilkj_l  Skl. 

(21) 

It  is  absolutely  clear  from  equation  (21)  that  each  iteration  is 
simply  now  a  series  of  matrix  operations  followed  by  a  matrix 
diagonal ization. 

The  next  step  is  to  perform  the  correlation  calculation. 
Consider  the  ordinary  second  order  Rayleigh-Schrodinger  case  here. 
The  second  order  correction  .to  the  energy  is  simply: 


E 


(2) 


=  2  2 
'i>j=o cc.  a>b=virt. 


(V 

i  jab 
_______ 

i  J 


-  V  >» 


i  jba 
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(23) 


The  coefficient  C  is  the  coefficient  of  the  basis  function  in  the 

P 

*th  Fock  orbital  for  the  p+n  configuration.  Thus  the  dominent 
correlation  problem  becomes  one  of  rotating  integrals  over  basis 
functions  to  integrals  over  Fock  orbitals.  This  again  is  simply  a 
series  of  matrix-like  steps.  General  considerations  for 
programming  equation  (21),  (22)  and  (23)  may  therefore  be  given. 

The  matrix  of  integrals  over  basis  functions,  gpqrs,  is 
sparse  (IX-  10X  density).  The  sparseness  has  two  causes,  one 
being  symmetry,  the  second  being  great  separation  of  basis 
functions.  Both  considerations  are  taken  into  account  before 
evaluating  any  integral,  thus  saving  time.  Once  the  integrals  are 
generated  one  need  be  more  particular  in  achieving  efficiency. 
Consider  as  an  example  the  part  of  the  Fock  matrix 

_(2)  _  _  n  —  P 

U  kJ  *ljkl  kl 

(22) 


or 


=  «XK  SK 


(22b) 


In  equation  (22b)  the  indices  i,  j  have  been  mapped  into  a  single 
index  I  and  k,  1  have  been  mapped  into  a  single  index  K. 

Computer  type  now  enters  into  our  consideration.  In  using 
vector  computers  like  array  processors,  Cyber  205' s,  Cray’s  or 


FPS164’s  one  should  avoid  logic  statements  inside  inner  do  loops 

in  most  cases.  Therefore,  the  sparseness  of  the  g  matrix  is  of 

little  help.  It  is  of  considerable  help  however,  on  conventional 

scalar — computers  and  one  may  work  directly  from  equation  (22)  on 

such,  loading  only  the  non-zero  integrals.  On  a  vector  machine 

however,  one  should  use  equation  <22b)  with  the  null  integrals 

included,  as  each  element  of  the  Fock  matrix  is  as  seen  a  simple 

vector  dot  product.  It  is  this  operation  which  is  maximally 

efficient  on  vector  computers  in  general.  Furthermore,  the  large 

2  4 

length  of  the  vector  in  (22b),  typically  of  length  10  to  10  ,  is 
ideal  for  such  systems  as  the  Cyber  205  as  well  as  the  other 
vector  machines.  Finally,  this  is  in  good  form  for  processing  on 
machines  with  parallel  architecture  such  as  a  Denelcore  or 
FPS164/MAX,  as  one  can  use  the  vector  Sk  as  a  constant  and  work  on 
several  of  the  Fj's  at  one  time.  Similar  considerations  apply  to 
the  matrix  operations  in  equation  (23).  Thus  code  developments 
are  machine  specific  to  ensure  maximum  efficiency. 

Recently  the  principal  investigator  has  used  the  method  of 
Kunz  and  Klein-*  to  incorporate  exact  boundary  conditions  at  the 
Quantum  cluster — Hades  interface  into  the  ICECAP  formalism.  This 
is  accomplished  as  follows.  Assume  in  the  region  beyond  the 
Quantum  cluster  boundary  that  the  first  order  density  matrix  for  a 
given  ion/atom/molecule  is  the  same  as  that  of  the  perfect  solid. 
This  first  order  density  matrix  may  be  trivially  determined  by  the 
method  of  local  orbitals.  A  program  to  evaluate  this  first  order 
density  matrix,  LOPAS,  has  for  a  reasonably  general  case  been 
developed  by  one  of  us  (ABK)  and  is  well  described  in  the 


literature.®  This  then  determines  the  entire  -first  order  density 
matrix  in  regions  outside  the  cluster.  The  entire  first  order 
density  matrix,  P,  is  formally  given  as 
*1  .  _  *  T  +1 


Ptxtx1  )  *  g  tx  x  )  +  tx^x  )+g  (x  x  ) 

not  in  c 


-  Pc|1  (^x1))  +meJ  P^,*1)- 

not  in  c 


(24) 


Here  P  refers  to  the  cluster  and  n  or  2  are  ion  sites  outside  the 


cluster.  Cross  terms  P^  ,  Pc^»  P^v  occur  because  in  Ref.  5  we 


exploit  the  generality  of  the  local  orbital  prescription  and  use 
different  localizing  potentials  for  different  regions.  The 
relevant  part  here  is 


P  <x  x1)  =  E  *.(x)*.<x1)S’.%  E  *  (x)*.  (x1) 
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Now  consider  the  general  local  orbitals  equation  for  defining  the 
*'s.  It  is  of  the  form 


CF  +  pwPl*i  = 


(25) 


F  is  the  Fock  operator  defined  in  Eq.  6  and  P  the  first  order 
density  matrix  defined  here  is  identical  to  that  in  Eq.  (5),  hence 


the  utility  of 
W  -  p 


this  method. 

V  S  . 

C  VE  pc’ 


For  in  c  we  arbitararily  chose 

(26) 


5 

Here  Pc  is  as  given  above  and  is  the  non  ionic  part  of  the 

potential  due  to  all  ions/atom/molecules  outside  of  the  cluster. 
Now  since: 


P  Pc  =  Pc  P  =  P 

one  finds 

'  'c  VE3  'c  "'i  ‘  £i  *i’  <27> 

becomes 

CF  ♦  Pc  VES  Pc>i  =  <28. 

This  is  exact  and  diefines  all  in  the  cluster  C.  Thus  the 

g 

added  term  Pc  VE  Pc  in  the  Fock  equation  replaces  the  crystal 

S 

outside  of  C.  Since  the  term  is  determined  in  this  method 

self  consistently  and  can  contain  even  the  contributions  of  an 
infinite  crystal,  this  approach  completely  solves  the  problem  of 
the  cluster  boundary  condition.  Current  efforts  are  directed  at 
implementing  this  equation  into  ICECAP.  Preliminary  tests  of 
ICECAP  on  moiecular/energetic  solids  should  be  made  early  in  1985. 


Some  Preliminary  Results 

F  and  Fft  centre  properties  (optical  absorption,  hyperfine 

constants,  and  spin  polarisation)  in  alkali  halides  were 

investigated  in  1980  by  Kung  et  al .  C13D,  with  a  nearest-nei ghbour 

UHF-SCF  cluster  based  on  free-atom  orbitals,  neglecting  lattice 
distortion  and  polarisation.  For  the  first  time,  the  Fft-centre 
absorption  splitting  was  described  correctly,  and  other 


encouraging  results  were  obtained  for  these  electrically  neutral 


defects,  although  the  smallness  of  the  cluster  was  recognized  as  a 
problem. 

Recently,  we  have  analysed  the  F+  centre  in  MgO  C141,  a 
charged  defect,  including  multipole  consistency,  based  on  HADES, 
and  using  ATMOL  C151  rather  than  Kunz’s  UHF  C4D  for  the  cluster. 

In  the  course  of  this  work  we  came  to  recognize  the  need  for 
complete-ion  pseudopotentials,  or  some  other  form  of  frozen 
extended  ions,  at  the  cluster  boundary.  Nevertheless,  the  crucial 
importance  of  correct,  consistent  treatment  of  polarization  in 
optical  absorption,  and  of  the  ion-size  effect  in  emission,  were 
demonstrated,  and  partial  electron  transfer  (hole  trapping)  in 
absorption  was  illustrated. 

These  methods  apply  directly  to  spectroscopic  computation.  A 
recent  series  of  such  studies  have  been  made  by  Goalwin  and  Kunz 
C163  for  atoms,  molecules,  and  solid  excitons.  In  the  case  of  the 
excitons,  the  change  in  electronic  structure  in  the  cluster  is  of 
a  quadrupole  nature  and  hence  outside  the  cluster  the  continuum 
dielectric  limit  is  adequate  and  included  here. 

The  results  of  this  investigation  are  tabulated  in  Table  i. 

The  energy  levels  calculated  in  the  UHF+MBF'T+polarization 
approximation  appear  to  be  within  .2  ev.  of  the  experimental 
levels  for  all  atomic  and  molecular  cases.  Singlet-triplet 
splittings  appear  to  be  accurate  to  within  .1  ev. 

One  point  frequently  brought  up  in  connection  with  the  UHF 
method  is  that  triplet  wavef uncti ons  are  not  eigenstates  of  Sz. 

Our  calculated  values  of  the  spin  in  the  triplet  state  were 
between  1.0000  and  1.0005  for  atomic  cases,  1.01  for  methane,  and 


between  1.600  and  1.002  -for  the  alkali  halides.  Since  Sz  is 
chosen  equal  to  1  it  seems  reasonable  to  identify  these  states  as 

triplet,  since  it  is  unlikely  that  a  mixture  of  eigenstates  of  S2 

with  substantial  contributions  from  states  with  eigenvalues  of  S= 

greatly  different  from  2  would  give  a  value  of  S  so  close  to  the 

expected  value  for  a  pure  triplet  state.  The  triplet  and 

contaminated  singlet  states  from  VHF  form  a  system  which  is 

redi agonal i zed  to  permit  formation  of  pure  singlet  systems  here. 

The  levels  hypothesized  to  be  the  result  of  a  3pe,4s  exciton  in 
the  alkali  halides  have  energies  and  singlet-triplet  splittings 
'dose  to  the  levels  theoretically  predicted  for  such  states.  The 
shoulder  in  the  NaCl  data  that  peaks  at  10.3  ev.  could  be  caused 
by  a  3pa3d  exciton,  as  the  theoret i cal  1 y  predicted  energies  for 
those  states  lie  in  the  region  covered  by  the  shoulder. 

The  splitting  of  the  calculated  energy  levels  for  the  3p°3d 
exciton  in  KC1  is  smaller  than  in  NaCl  and  the  observed  peak  in 
KC1  is  narrower.  However,  the  energy  observed  is  .7  ev  to  1.0  ev. 
higher  than  the  calculated  energy.  It  is  possible  that  this  is 
due  to  the  unusually  low  correlation  energy  of  about  0.04  Hartrees 
obtained  with  this  basis  set.  It  would  be  difficult  to  increase 
the  size  of  the  basis  set  and  still  do  the  calculation  on  the 
computers  available  to  us.  It  is  also  possible  that  the  computed 
3d  level  refers  to  the  first  d  level  possible  in  the  excited  state 
as  the  present  use  of  the  variational  principle  would  imply.  In 
such  a  case,  the  narrow  peak  identified  as  a  d-level  in  the 
absorption  spectra  could  well  be  an  accumulation  point  for  the  p 
to  d  transition  series  limit.  It  seems  to  us  that  the  latter 
explanation  is  likely. 


v.v.v. 


No  variational  collapse  -for  the  excited  states  is  observed  in 
any  of  these  systems.  Our  computed  value  for  the  He  ls2s  singlet 
energy  is  0.7562  Hy.  Fraga  and  Briss  C17D  obtain  .7461  Hy.  Other 
workers  C18,  191  have  not  reported  the  ground  state  energy  they 
obtain  but  give  absolute  energies  for  the  excited  state  comparable 
to  our  result. 

Polarisation  and  Distortion  of  the  Surrounding  Lattice:  The  HADES 
Code 

The  HADES  code  is  based  on  the  general  formulation  for 
treating  the  defective  lattice  developed  by  Lidiard  and  Norgett 
(20)  and  Norgett  (21).  It  relies  on  the  idea  that  the  total 
energy  of  the  system  is  minimized  by  a  relaxation  of  the  ions 
surrounding  the  defect  and  that  this  relaxation  decreases  rapidly 
as  the  distance  from  the  defect  increases.  We  may  therefore 
partition  the  crystal  into  an  inner  region  I  in  which  the  lattice 
configuration  is  evaluated  explicitly  and  an  outer  region  II  that 
can  be  viewed  from  the  defect  as  a  continuum.  The  total  energy  of 
the  system  may  then  be  written  as 

E  =  EI(^  +  Eint(^>  +  EII<*> 

in  which  Ej  (x.)  is  the  energy  of  the  inner  region,  Ejj.(f)  the 
energy  of  the  outer  region  and  ^  <x_,  {.)  the  interaction  energy 
between  regions  I  and  II.  x.  is  a  vector  of  the  independent 
coordinates  describing  the  inner  region  and  £  is  a  corresponding 
vector  of  the  displacements  in  the  outer  region.  £  is  both 


formally  distinguished  from  x.  and  assumed  to  be  an  implicit 


approximation. 

We  may  eliminate  the  outer  crystal  term  Ejj<$)  by  assuming 
that  the  outer  region  responds  harmonically  and  so  Ejj(J)  is  a 
quadratic  function  of  Applying  the  equilibrium  condition 


then  allows  us  to  eliminate  Ejj(f). 

In  principle,  we  may  now  find  the  defect  energy  by  applying  a 
force  balance  condition  to  the  ions  in  the  inner  region  i.e. 

HI  -«>' 

t  =  constant 

To  proceed  further  we  require  an  explicit  representation  of 
the  energy  E.  We  assume  this  to  be  the  sum  of  two-body 
interactions.  Details  of  how  this  is  implemented  in  the  HADES 
program  and  of  the  minimization  procedures  used  to  obtain  the 
force-balance  position  are  given  in  references  (21)  'and  <22) . 

The  HADES  code  and  its  extensions  have  been  used  as  the  basis 
of  several  hundred  papers  concerning  defect  behaviour  in  ionic 
crystals.  Extensions  to  treat  defects  at  surfaces  and  interfaces 
(like  grain  boundaries,  etc.)  already  exist,  and  may  be  used  in 
later  phases. 

Some  specific  progress  in  implementing  these  ideas  for  systems 
such  as  CHa,  (NO) a  and  nitromethane  has  been  achieved  and  will  be 
reported  in  an  interim  technical  report  shortly.  Nonetheless,  the 
general  thrust  of  the  past  year  has  been  directed  toward  a  general 
implementation  of  these  rather  comprehensive  ideas  for  studies  on 
defects  and  impurities. 
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1.  Energy  Trapping,  Release  and  Transport  in  Three  Dimensional 
Energetic  Solids,  A,  Barry  Kunz,  Principal  Investigator, 
Office  of  Naval  Research,  awarded,  N00014-81-K-0620,  1  July 
1981-30  June  1986,  total  award  $345,768.00.  Current  annual 
level ,  $75, 000. 00. 

2.  F'ropert i  es,  of  Bound  States  of  Transition  Metal  Negative  Ions 
D.  R.  Beck,  Principal  Inverstigator ,  National  Science 
Foundation,  pending,  1  October  1985-30  September  1987,  total 
request  $58,041. 

3.  Des'elopment  of  Ab-Initio  Molecular  Potentials  for  Certain 
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Gas  Research  Institute,  pending,  1  January  1986-31  December 
1988,  total  request  $142,500. 
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